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We investigate the inelastic hard disk gas sheared by two parallel bumpy walls (Couette-flow). In 
our molecular dynamic simulations we found a sensitivity to the asymmetries of the initial condition 
of the particle places and velocities and an asymmetric stationary state, where the deviation from 
(anti)symmetric hydrodynamic fields is stronger as the normal restitution coefficient decreases. For 
the better understanding of this sensitivity we carried out a linear stability analysis of the former 
kinetic theoretical solution [Jenkins and Richman: J. Fluid. Mech. 171 (1986)] and found it to be 
\ unstable. The effect of this asymmetry on the self-diffusion coefficient is also discussed. 



in 



o 
o 



> 



c3 

B 
-6 



45.70M,45.70Q,51. 10,51.20 



O ' I. INTRODUCTION 

In the last decade the clustering instability of undriven granular gases was extensively researched. Since the first 
explanation of the instability Jffll investigation had been done with several methods including stability analysis of 
the hydrodynamic equations [gfl3[] and fluctuating hydrodynamics . The examination of the influence of nonlinear 
coupling between hydrodynamic modes was also performed j^], and the long time behavior of the clustered state 
was with mode-coupling theories examined ||. The evolution of vortex velocity patterns preceding the clustering 
instability is well understood 0, and the minimal system size, where the instability appears, and its dependence on 
the restitution coefficient are revealed . 

A driven configuration, the uniformly sheared inelastic hard sphere gas (with periodic boundary conditions) shows 
also an instability and a resulting pattern formation. In S it was shown that this pattern formation was guided by 
an instability for the short times and a shearing caused convection for longer time-scales. The special importance of 
this system is the coupling of macroscopic and microscopic length-scales in sheared situations as discussed in P-pH . 

A more realistic configuration is the Couette-flow, the shearing of the gas between two parallel walls moved in 
opposite directions. For this configuration with walls consisting of disks Jenkins and Richman derived boundary 
conditions for the momentum and heat transfer of the walls [O] proceeding from their kinetic theory for dense granular 
gas es |l3| ] . In subsequent years this theory was developed further p"4| ] for non Maxwellian velocity distributions as in 
jp^l . In all of these papers the resulting hydrodynamic fields are symmetric (density and granular temperature) or 
antisymmetric (flow field) on the half-point between the two bounding walls. In view of the clustering instability it 
is of major interest whether this instability occurs in this configuration and, if so, what are its consequences in the 
Couette-flow. This paper is devoted to the investigation of these questions. 

The paper is organized as follows: In the next Section (|J) we present molecular dynamic simulations of this system 
where we found an interesting sensitivity to the initial condition of the simulation and asymmetric hydrodynamic 
fields contradicting the results in JOs]. Therefore we carried out the linear stability analysis (Sec. |l|) around the 
O ■ solution of Jenkins and Richman |12| and found it to be unstable against certain fluctuations. We discuss the effect 
of this instability on the diffusion coefficient in Sec. IV and close the paper with summary in Sec. 

•i-H . 



II. SIMULATION RESULTS 

The system considered here consisted of N identical, inelastic hard disks with mass m = 1 and radius r = 1 confined 
in a rectangular area of the size L x x L y . This two-dimensional area is bounded by two parallel walls from two sides 
which also define the direction x and are of the length L x . The system is closed through periodic boundary conditions 
in the x direction. The walls are L y distance apart (y direction) . The origin of the coordinate system is placed in the 
middle of the simulated system what sets the center of the wall disks to y — ±L y /2. 

The disks interact through the inelastic rough hard sphere potential meaning instantaneous two particle collisions 
characterized by given ratios of the final and initial velocities in the normal and tangential direction as the normal 
and tangential restitution coefficients e„ = —v^/v^ and e t = v{ jv\, f and i meaning the final and initial velocities. 



1 



In case of sliding contacts et is replaced by the ratio of tangential and normal momentum transfer characterized by 
the friction coefficient \i. As a driving force we move one of the walls with constant velocity 2U. 

We investigated the system with the event-driven molecular dynamics method ideal for simulating instantaneous 
collision rules jig] . Characterizing the collisions we used et — —0.3 and /i = 0.2 and varied e n between 0.6 and 
0.99. We consider the system sized L x = L y = 40 and the wall velocity 2U — 10 if not otherwise mentioned. The 
post-collision velocities and angular velocities as a function of the pre-collision ones are (for m — 1, r = 1), 

fj)r n )r n + — r t (1) 
Vj)r n )r n + — r t (2) 
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r n being the normal vector of disk surfaces in collision r n — (rj — Ti)/\rj — r, | and r t the tangential vector, the normal 
vector rotated with 7r/2 in anti-clockwise direction. Ii, Ij are the moments of inertia of the disks. The parameter j t 
depends on the type of the collision in tangential direction, 

1 + e 

jt = -/x — 2~~^( v ' ~ V ^ r " if ^ > ^° ( slidin g contact); (5) 
jt — 6t - Vt if (f> < (f>o (sticking contact); (6) 

where v t is the tangential velocity 

v t = (v; - Vj)r t + r(uj + utj) (7) 
and 0o separates the sliding and sticking regions in the angle of incidence <f> — |vt|/((vj — Vj)r„), 

= ■ (8) 

2 1 — e* 

To avoid inelastic collapse |l6f| , the intrinsic numerical breakdown of the method, we stopped the simulations 
when the time interval between two subsequent collisions became smaller than the precision of the computations as 
proposed in Q|; but this stop occurred only for the used smallest restitution coefficient e„ = 0.6 and for special initial 
conditions where most of the particles did not move, and the well-known chain-like arrangement of particles could 
evolve after the simulation started. For the experiments below we used particle numbers N = 100 — 240 meaning 
area fraction v £ [0.245,049]. As test runs for larger systems we used {L x = L y = 100, JV = 500, {y — 0.162)} and 
{L x = L y = 100, N = 1000, (v — 0.324)} parameter sets resulting in the same phenomenon. 

Using several different starting configurations we noticed that there appears a dense region at one of the walls in 
the system which is more and more distinct with decreasing e n and not recognizable for e„ near 1, FIG 1. Using 
randomly placed particles and an initial uniform velocity distribution with several mean values in the x direction we 
observed that the initial value of the mean particle velocity determines the wall which is chosen for building up a 
clustered regime. If we define Vd as 

v d = (v) - U 

we can characterize the final place as if Vd > the upper y = L y /2 wall is chosen, if v d < the lower y = —L y /2 one, 
FIG 2. We used different initial conditions to test this finding jTif . With particles placed in a stripe in the middle of 
the system, organized on a lattice with initial velocity U , and using two of them shot against the two wall with initial 
velocities which were Galilei-symmetric with the situation of two walls moved with ±Z7 we could maximize the time 
needed to develop the non-symmetric density field but it appeared also after this initial condition because numerical 
errors provide the needed fluctuations to drive the system in one of the steady states. 

We did simulations for smooth disks (et = 1 and [i = 0) as well and found the same effect. Namely, if we chose 
e„ to be small enough there is a distinct band of particles at one of the walls and for e n close to 1 this band does 
not evolve. The normal restitution below which we can speak about an asymmetric phase depends on the density 
as for higher densities it appears at higher e„ values on the same system size L y . It also depends on L y for fixed 
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density as for larger L y the kinetic energy influx per particle decreases as it is proportionate with the length of the 
wall. The assumed critical restitution coefficient e c (p, L y ) does not depend on the driving velocity as U is the only 
parameter which includes dimensionality of time and changing U is a rescaling of the time unit. For the investigation 
of intermediate e„ values for which the dense stripe is not so definite we measured density and velocity profiles and 
the granular temperature in the system and found the same effect now in the asymmetric nature of this functions 
FIG 3. One can also observe that with decreasing restitution the minimum of the granular temperature moves to one 
of the walls. 

We found these results valid for a broad range of densities and for different system sizes. For low restitution there 
appeared a crystalline structure as compact layers near a wall. The fact that the simulations did not stop at such 
dense configuration implies that the disks moved around a fixed place and did not form a long enough connected line 
of particles in which inelastic collapse could occur [|l6| . 

Simulations with time step driven molecular dynamics using damped harmonic oscillator forces between the particles 
as a function of their overlap 1 15| - for achieving velocity independent restitution - show the same clustering effect at 
the walls. This suggests that this effect is closely related to the dissipative collisions and does not require instantaneous 
collision rules. 



III. STABILITY OF THE JENKINS-RICHMAN SOLUTION 

Our above findings contradict the kinetic theoretical calculation of Jenkins and Richman ]T^| and the improved 
versions of the problem |jl4| inasmuch as the cited results are always (anti)symmetrical profiles in the y direction for 
the hydrodynamic fields of density, flow velocity and granular temperature. This contradiction raises the question of 
the stability of the (anti)symmetric solution. Moreover, it is plausible that in an elastic system the hydrodynamic 
profiles are symmetric, though this is not a steady state because of the ever rising temperature caused by viscous 
heating. Therefore it is of interest to see weather a phase transition occurs in the system, specifically, weather there 
exists a critical e n value or the (anti)symmetric solution becomes unstable at arbitrarily small inelasticity. Because at 
high densities the instability in question arises at normal restitution closer to 1 and kinetic theoretical calculations are 
more and more inadequate for increasing inelasticity it is reasonable to use kinetic theories for dense systems and we 
desist from the consideration of more sophisticated hydrodynamic equations for low density systems |l8],[l9| . Another 
reason for taking theories for dense gases into consideration is that the underlying assumptions are more valid as the 
ratio of the simulated system size to the mean free path is much larger in a dense gas. For this reasons we performed 
a linear stability analysis of the Jenkins-Richman solution. We now first recite the equations of Jenkins and Richman 
p3[ and the boundary conditions and the solution of the problem ]T^ ]. Then we perform a linearization around this 
solution and calculate numerically the eigenvalues of the stability matrix. 



A. Jenkins-Richman Solution 

In this section we briefly describe the Jenkins and Richman p2| solution for a two-dimensional system of inelastic 
hard disks driven by two parallel bumpy walls. The hydrodynamic equations for the density p, flow velocity u and 
fluctuation energy T for smooth disks without external forces are according to |l3| : 

P = -pV-u, (9) 
pu = -V • P , (10) 
pf = -V • Q - Tr(P ■ Vu) - 7 . (11) 

Here the dissipation rate 7 and heat transport coefficient n are of the form, 

7 = 4(i -:- )Kr (12) 

2pcrvg Ti 

K= T (13) 

7T2 

where go is the Enskog correction term accounting for excluded volume effects as a function of the packing fraction v: 

16 -7u 
30 ~ 16(1 -v) 2 ' 
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and a = 2r is the disk diameter. The constitutive relations for the energy flux Q and pressure tensor P are, 



Q = -kWT (14) 
P = (2pvg T - i«Tr(D)I - kD. (15) 

From the boundary geometry and the assumption of a Maxwellian velocity distribution the momentum and fluctuation 
energy supplied by the wall in unit time and length are determinable This needs an expansion in e = a jL y and 
assumptions on the dependence of given ratios of the hydrodynamic fields on e (see [p"2|), what narrows the validity 
of the theory. The supplied momentum obtained that way is 

M a = \p X {l + e n )T L a +(IY (-^ - cosfl) + 0) ' ^r^I a ^) (16) 



with 



(2 \ 2 

Ia/3-y = ( g sin 9 — 2 J n Q nan 7 — — sin 9 (riatpty + n^t a t 1 + n^tatp) (17) 

and the supplied fluctuation energy is 

D= (2_\i px (l-e)Ti 
\7ry sin 9 

where n and t are the normal and tangential vectors of the wall respectively and 9 characterizes the bumpyness of 
the wall. For wall disks with the same disk diameter as of the gas particles it equals |l^] : 

sin 9 = — 

a 

where d is the distance separating the centers of two wall disks. For the wall geometry of our considered system 
d — 2a and therefore sin 9 = 1/2 and 9 = tt/6. In ( p^ , p^| ) \ accounts for static correlation effects in collisions at the 
boundary and v a , a = x, y is the slip velocity, the difference between the velocity of the wall and the flow v = U — u. 
From the equality of momentum and energy transfer of the wall for unit time and unit length the boundary conditions 
are: 

M = P-n, (19) 
M • v D = Q • n . (20) 

The solution can be calculated from these equations except one parameter - the value of solid fraction at the boundaries 
(for fixed number of particles) - which must be iterated to get the given density in the system what we did numerically. 
Thus the form of the hydrodynamic fields solving the boundary value problem is, with T(y) = w(y) 2 , 

X(U-v)Na , f\y\ 
w(v) = Z i . —-77——— cosh I -p ) , (21) 



2tt2 smh(X/2)SL \ L 
«V) = ^K<**(¥) (22) 



(23) 



sinh(A/2) V L 

with S and N being the shear stress and pressure respectively and A must be solved from [JL2| 

,2 



A tanh /A\ _ (1 - e n )L9 I sin0[l - (4V2^/3cr) sin 2 9} 
-tan - I _____ 



The density, calculated from 



X = 2^5o (24) 



must fulfill the condition that its integral on the whole system gives the prepared particle number. As mentioned, in 
that condition the only free parameter remains the density at the boundaries what must be calculated iteratively. 
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B. Stability Analysis 



Around the above solution of the boundary value problem we performed a linear stability analysis with the perturbed 
quantities, 

pi = p + 8p (25) 

u' = u + Su (26) 

T' = T + 5T. (27) 

We split the velocity perturbation in the x and y direction and use plane waves for the perturbations, 

5p = 5pkexp(~ik ■ r) (28) 

' U Wf k W-&.r) (29) 

ST = 5T k exp(-ik ■ r). (30) 

With this choice the resulting equations have the form 



/ d t Sp k 

d t 5w k 
V d t ST k 




(31) 



where A is the stability matrix with the elements given in Appendix |a|. 

With the help of MAPLE V Release 5 we analyzed the eigenvalues of the stability matrix A for different densities, 
restitution coefficients, wave numbers and sites. On FIG 4-7 we present results for some of the parameter sets which 
are representative for the different situations. With the choice of the negative sign in (^l|) a negative eigenvalue means 
an unstable solution. All of the shown figures feature at least one eigenvalue which is negative in the small wave 
number region. This signalizes the unstable fluctuation. The upper zero point of this eigenvalue nears to the zero 
wave number as we increase the restitution coefficient and becomes smaller than 2tt/L 1 the smallest possible wave 
number in the system, at a certain value of e n what depends on the other parameters. However, this value is very 
close to 1 what suggests that the symmetrical solution is always unstable linearly and nonlinear effects can move this 
transition point to higher inelasticities. We analyzed the most unstable direction of the solution also as we plotted the 
unstable eigenvalue for a given amplitude of the wave number and for the whole angle measured from the direction 
of the mean flow in anti-clockwise direction as seen in FIG 8,9. According to the figures the most unstable directions 
are close to the ±7r/4 angles which are the same angles found to be the governing unstable directions in linear order 
for the pattern formation by Tan and Goldhirsch for the uniformly sheared granular gas M. 



IV. SELF-DIFFUSION COEFFICIENT AND CLUSTERED FLOWS 



In this section we investigate the crossover to asymmetric stationary states of the system and try to show the 
effect of the instability discussed above on the transport coefficients. For this reason we measured the mean square 
displacement of the particles, 

rl = ((x(t) - (v x )t - x(0)) 2 ) , rl = {(y(t) - y(0)) 2 ) (32) 

and velocity correlation functions, 

c x = ((v x (t) - (v x ))(v x (0) - (v x ))) , cy = (v y (t)v v (0)) (33) 

in the system for stream-wise and perpendicular directions. In ( |32"| ) and (^) the brackets () mean space and ensemble 
averages over the whole system. (v x ) was calculated from the measured average displacement of the particles averaged 
over several (50-100) runs and over particles. Ensemble averages were carried out the following way. For every 
parameter set we let relax the system in its steady state. Then measured granular temperature averaged over stripes 
of width 1 = r/2 generating the T(y) function. Before every measurement we perturbed the velocities of every particle 
with a uniformly distributed velocity the support of the distribution being 2y/T{y), where y was the coordinate of 
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the particle. After relaxation we carried out the measurement. The (v x ) average velocity was evaluated from a linear 
fit on the averaged x displacement of the particles being the steepness of x as a function of time. 

From the mean square displacement for long times we could evaluate the diffusion coefficient of the particles. We 
fitted a power law with two fit parameters D and b, 

/(*) = 2Dt h 

on the mean square displacement restricted to long times. If the value of b was up to errors (±0.02) near 1 the 
system was assumed to be in a diffusive state. After that we took 6=1 and fitted f(t) again for D. The time the 
system evolves into this diffusive state depends on e n as suggested in [|o| for the homogeneous cooling state, and also 
depends on the steady state configuration which is also e„ dependent. This has several reasons. First, to achieve 
this diffusive behavior in x direction the particles have to experience the characteristics of the system - finiteness 
and inhomogeneity of the hydrodynamic profiles - therefore have to travel between the to walls several times. As e ra 
decreases and a dense cluster evolves at one wall this recognition takes more and more time. Also the form of the 
velocity correlation function shows increasing correlations with decreasing e n as expected. After a short (t w 10) fast 
decay c x shows an exponential decay up to times t ps 100 and changes to a slow decay for later times (small e„) or is 
relaxed (large e„) FIG 10. Both the short-time correlations and the characteristic time of the exponential increases 
as we decrease e n , along with the earlier changeover to the slow decay if it is apparent. 

The plot of the mean square displacement versus time diagram (FIG 11) shows a plateau value of the diffusion 
coefficient for larger values of e„. Below e c = 0.82 D increases fast with decreasing e„. We suggest that this value 
marks the transition point where the particles begin to prefer one of the walls and the appearing dilute regions allow 
for a greater mobility for the particles. Making use of the noticed exponential decay of the velocity autocorrelation 
c x we also estimated D as the integral of a fitted exponential, 

D w At if c x fa Ae~ t/T . 

knowing that it underestimates D as Aer 1 ^ is smaller than c x for very short and for long times, but it gives nearly the 
same result as the long time behavior of r 2 x above the transition point and shows a smaller increase below it (FIG 11). 
This breakup comes from the fact that under the transition point there appears a slow decay in the velocity correlation 
function after the exponential decay used to approximate it. The short-time exponential decay characterizes only the 
dilute part of the system. The slow decay thereafter comes from the averaged effects of the exigous diffusion in the 
dense region and from the influence of particles entrapped or emitted by the dense region from or to the dilute one. 
In absence of the asymmetric dense state above the transition point e c only the exponential relaxation appears. 

The buildup of correlations is also a cause for increasing times needed to reach the diffusive state thus increasing 
simulation time what restricted us to smaller system sizes for dense material as it did not allow fast computations. 
Therefore the investigation of the size and density dependence of the transition point is left for later studies. 

The retrievable information from the mean square displacement in y direction is bounded as the system itself namely 
it goes to a plateau value according to L y . This finiteness influences velocity correlations also. c y shows only short- 
time fast decaying correlations (and there appear anti-correlations for small e n .) With the methods discussed above 
we do not noticed any mark of the transition (only a small sign as presented in pi[ for time-step driven simulations). 

For more definite insights in the behavior of the system we carried out measurements specific to the configuration 
of the stationary states. We divided the system into stripes parallel to the walls and restricted the averages in ( |33| ) 
on particles starting from these stripes, 

4 = <(« a (t) - vi(t)) {v x (0) - 4(0))), , 4 = (v y (t)v y (0)) i (34) 

where <>i denotes averaging over particles being in stripe i at t — and the ensemble average described above. The 
velocity v l x (t) is the average velocity of particles starting from stripe i at time t. This velocity correlation function 
characterizes the stripe more and more as the difficulty to leave the stripe increases which is the case in the dense 
layer evolving at the wall. The fact that the velocity correlation function c x does not show a negative minimum at 
short times (FIG. 10) - which would be a sign of caging effects - suggests, that the particles essentially comoving 
with the wall move barely in the x direction compared to one another. Instead because of the sheared situation they 
are driven by diagonal collisions which dominate the motion and smooth out c x . There is a negative minimum in c y 
(0 indexing the stripe nearest to the wall at the cluster) but this includes the effect of the collisions with the wall 
particles which give a negative contribution as also seen on the velocity correlations in the stripe at the opposite wall. 
This appears also at restitution coefficients well above the transition. 
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V. DISCUSSION 



In this paper we considered a granular gas sheared by two bumpy walls. By means of molecular dynamics simulations 
we discovered that the system shows a spontaneous symmetry breaking as it evolves to a stationary state which 
is asymmetric in the hydrodynamic fields with respect to the centerline of the system. The location where this 
density maximum (temperature minimum) is shifted to from the symmetry axes is predetermined by the initial 
velocity distribution as the system chooses the wall according to this distribution as discussed in Sec. 0. During the 
simulations, measuring the mean y coordinate of the particles, we did not observe any switch between the possible 
two states where the asymmetry was explicit. We found in a linear stability analysis that the symmetrical solution 
proposed by Jenkins and Richman [jl2| is unstable. In analyzing the solution one finds that the description of the 
asymmetric profiles would need higher order equations in e as in . We do not consider this analysis as a strict result 
quantitatively. According to p2[ the equations of Jenkins and Richman are not suitable for linear stability analysis 
because they neglect the time dependence of the mean free path. But the equations with the corrected constitutive 
relations jig ] have the same dependence on the hydrodynamic field and differ in the e„ dependence only. We assume 
that this correction terms cannot stabilize the system but can at most increase the e„ value under which the solution 
becomes unstable. 

We found a connection to the unstable directions observed in uniformly sheared granular gases j8[. That suggests 
that not uniformly sheared systems possess the same instability but in the Couette-case the pattern formation is 
repressed by the boundary conditions. However, the geometry of the boundary conditions (walls) plays a major role 
in stabilizing the flow fields leading to a stationary state. In the boundary conditions of Jenkins and Richman Jl^] 
only static correlations are considered and dynamic correlations, the role of multiple and correlated re- collisions with 
the wall, are neglected. We measured the number of collision sequences which consist of collisions of a disk with 
two wall disks successively because such collisions provide disks leaving the wall in a steeper angle what can provide 
greater contribution to the pressure at lower granular temperature. In the asymmetric flow regime the number of such 
collisions is clearly higher at the preferred wall than at the other what supports our assumption. This is also implied 
by the fact, that in the final stationary state the minimum of the granular temperature is only at the wall for very high 
inelasticities or densities where we could observe closed layers of particles at the wall chosen. For moderate densities 
and restitution the temperature minimum remains in the bulk of the system, and the clustering process stops at a 
stable density profile. Such effect was described for a constrained homogeneous system in ||. This boundary effects 
are one reason for the different results for the transition point obtained from the stability analysis or the measurement 
of the self-diffusion coefficient and velocity autocorrelation functions respectively. The found effect in the dependency 
of these functions on the normal restitution coefficient described in Sec. [[^ is well described as a consequence of the 
transition. However, the system size used in the simulations is too small to compare the two results. We suggest that 
for larger systems where the Jenkins-Richman theory is more valid we would obtain larger values for the transition 
point but the involved measurements leave this for later studies. As a final note we would like to mention that in 
systems with elastic collisions with the wall disks the flow with the same parameters remains symmetric. 

The consequences of this paper are hard to verify in experiments because of the absence of gravity. However 
microgravity experiments for the same configuration are under way ^3|. We suggest that the two-dimensional case 
considered here would be realizable with disks on an air table as this would exclude gravity effects and would not 
interfere with the behavior of the disks in the crucial horizontal direction. Periodic boundary conditions could be 
achieved with two similar stadion-like chains of disks as performed in |23(] . 
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In the following we describe the matrix elements of the linear stability matrix. The notations include p = "Qv where 
•d = 7r/4 for a = 1, m = l and ip — d y (vgQ(v)). 
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FIGURE CAPTIONS 

FIGURE 1: Snapshot configurations in the steady state for different restitution coefficients: (la) e„ = 0.7, (lb) 
e„ = 0.8, (lc)e„ = 0.9 Lines from centers of particles indicate the direction and magnitude of its velocity. 
FIGURE 2: Snapshot configurations in the steady state for different initial velocity distributions (e„ = 0.7). 2a: 
v d < 0, 2b: v d > 0. 

FIGURE 3: Granular temperature in the x and y direction as a function of y (in units of disk radius r = 1) across 
the system for e„ = 0.9. Values are averages over stripes of width 4 parallel to the walls with diameter 4 for a system 
of size L y = 40. 

FIGURE 4: Eigenvalues of the stability matrix for e„ = 0.9, y — 0, average packing fraction v — 0.486, and k a points 
in the direction a = ir/4 measured anti-clockwise from the mean flow direction. 

FIGURE 5: Eigenvalues of the stability matrix for e„ = 0.95, y = 0, k x = 2-k/L and average packing fraction v — 0.486 
as a function of k y L/ir. 

FIGURE 6: Eigenvalues of the stability matrix for e n = 0.997, y = 0, average packing fraction v = 0.486, and k Q 
points in the direction a — ir/4 measured anti-clockwise from the mean flow direction. 

FIGURE 7: Eigenvalues of the stability matrix for e„ = 0.999, y — 0, average packing fraction v = 0.486, and k Q 
points in the direction a = 7r/4 measured anti-clockwise from the mean flow direction. 

FIGURE 8: Angle dependence of the real part of the eigenvalues for e n — 0.95, y = and v — 0.486. The eigenvalue 
is plotted as a function of the angle measured from the flow direction k x = (2n/L) cos(0) and k y = (2ir/L) sin(</>). 
FIGURE 9: Angle dependence of the real part of the eigenvalue for e„ = 0.997, y — and v — 0.486. The eigenvalue 
is plotted as a function of the angle measured from the flow direction k x = (2ir/L) cos(</>) and k y = (2n/L) sin(</>). 
FIGURE 10: Velocity autocorrelation functions c x (in units of U 2 = 25) for a system with parameters L x = 20, 
L y = 20, N = 50, {v = 0.426). In order of increasing values for short times are e„ = 0.8,0.78,0.77,0.76. Inset: 
Velocity autocorrelation functions c x (same units) for a system L x = 80, L y = 20, N = 200, [y — 0.426) and e = 0.76. 
The system was divided in 5 stripes of width 4. Correlation functions are shown only in the stripes at the walls. 
Solid line marks the function at the wall of the cluster, and dot-dashed line at the opposite wall. Time is measured 
in natural units (r = 1, U = 5). 

FIGURE 11: The self-diffusion coefficient as a function of the restitution coefficient e„ for L x = L y — 20, N = 50, 
(D = 0.426) measured in natural units. This runs were carried out moving both walls in opposite direction with 
U = ±5. Diamonds (O) marking values determined from the mean square displacement, and triangles (A) marking 
values determined from the velocity autocorrelation function. 
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